A two-dimensional paradigm for symmetry breaking: 
the nonlinear Schrodinger equation with a four-well potential 
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We study the existence and stability of localized states in the two-dimensional (2D) nonlinear 
Schrodinger (NLS)/Gross-Pitaevskii equation with a symmetric four- well potential. Using a four- 
mode approximation, we are able to trace the parametric evolution of the trapped stationary modes, 
starting from the corresponding linear limits, and thus derive the complete bifurcation diagram for 
the families of these stationary modes. The predictions based on the four-mode decomposition 
are found to be in good agreement with the numerical results obtained from the NLS equation. 
Actually, the stability properties coincide with those suggested by the corresponding discrete model 
in the large- amplitude limit. The dynamics of the unstable modes is explored by means of direct 
simulations. Finally, while we present the full analysis for the case of the focusing nonlinearity, the 
bifurcation diagram for the defocusing case is briefly considered too. 

I. INTRODUCTION 



In the recent years, there has been a considerable effort on experimental and theoretical studies of Bose-Einstein 
\ condensates (BECs) pj, [2|. Many of these studies were focused on macroscopic nonlinear structures that arise in 
p I . BECs, which often have counterparts in nonlinear optics [3]. One of the appealing features of this setting is the 
possibility to tailor the desirable geometry of magnetic, optical, or combined traps that confine the ultracold bosonic 
atoms. For this reason, the analysis of the existence, stability and dynamical properties of nonlinear localized states 
trapped in these geometries have become a focal point of research. The theoretical analysis is enabled by the fact 
that a very accurate description of dilute atomic BECs is furnished, in the mean-field approximation, by the Gross- 
Pitaevskii (GP) equation, which is a variant of the nonlinear Schrodinger (NLS) equation. The cubic nonlinearity in 
the GP equation originates from the interatomic interactions, accounted for through an effective mean-field. The NLS 
equation in this, as well as in somewhat different forms, is relevant to a variety of alternative physical applications in 
nonlinear optics and other areas [3, 0, [5[ . 
(S| ■ Among the trapping configurations available in current BEC experiments, one that has drawn considerable attention 
' is the double- well potential (DWP). Its prototypical realization is provided by the a strong parabolic (harmonic) trap 
\ combined with a periodic potential, which can be created as an "optical lattice", by a set of coherent laser beams 

■ illuminating the condensate P, H, The DWP was realized experimentally in [6], using the magnetic field to 
\ induce the parabolic trap. The experiments reported in Ref. [6] revealed a variety of fundamental effects, including 

■ tunneling and Josephson oscillations for a small number of atoms, or macroscopic quantum self-trapping leading 
to a stable asymmetric partition of atoms between the wells for a sufficiently large number of atoms. DWPs have 

^. J . also inspired theoretical studies of various topics, such as finite-mode reductions, finding exact analytical results for 
/\ ■ specially designed shapes of the potential, quantum effects [7, 8, 0, [lO^Iiil^Illl^E^H^IlBI, and a nonlinear DWP (alias 
double- well pseudopotential), which is induced by the respective spatial modulation of the nonlinearity coefficient [16|. 
It is relevant to mention that DWPs have also been studied in the context of nonlinear optics, including twin-core 
self-guided laser beams in Kerr media [17] and optically-induced waveguiding structures in photorefractive crystals 

The aim of the present work is to extend the analysis of the DWP to a two-dimensional (2D) setting. Unlike previous 
studies of the trapping of quasi-2D BECs under the combined action of harmonic traps and optical-lattice potentials 
[195, we implement a Galerkin-type few- mode reduction to deduce a discrete model, based on a symmetric set of four 
wells, which is the most natural configuration in the 2D case. The same model can also be realized in optics, using a 
bulk nonlinear medium with a set of four embedded waveguiding rods. In that setting, we use numerical methods to 
generate a bifurcation diagram for possible stationary states of the system. It is worth noting that all the states that 
are expected on the basis of a four-site discrete nonlinear Schrodinger (DNLS) reduction [20] are also obtained in the 
continuum model with the combined parabolic and periodic potential considered herein. Furthermore, their stability, 
in the large nonlinearity limit, coincides with what is expected from the DNLS model. 

The paper is structured as follows. In section [III we present the model and the derivation of the four- mode 
approximation. Numerical results are reported in section lllli We present complete bifurcation diagrams of the 
possible stationary states for both the underlying GP equation and for its four-mode reduction. Comparison between 
them demonstrates very good agreement. In addition to the study of the existence and stability, evolution of unstable 
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modes is explored too by means of direct numerical simulations. Finally, we summarize our findings in section HV} 
where we also discuss possible directions for further work. 

II. THE MODEL AND THE GALERKIN APPROXIMATION 

We start by presenting the basic model in the quasi-2D setting, namely the NLS/GP equation in (2 + l)-dimensions, 
which is expressed in the following dimensionless form [1, ^, ^] : 

idfU = Lu -\- s\ufu — jiiu, (1) 

where u{x, t) is the normalized wave function, ji the chemical potential (— /i is the propagation constant in the optical 
realization of the model), 5 = ±1 corresponds, respectively, to repulsive or attractive interatomic interactions in BEG 
(alias self-defocusing or self- focusing Kerr nonlinearity, in terms of nonlinear optics), and L is the single-particle 
operator given by: 

L=-^A + T/(x,y). (2) 

In Eq. ([2j), A = 9^ + 9^ is the 2D Laplacian, while V{x^y) is the four- well potential, assumed to take the following 
form, 

V{x, y) = il^V^ + Vo [cos{2kx) + cos{2ky)] , (3) 

with = x'^ -\-y'^. It is clear that V{x^ y) is composed of a harmonic trap of strength and a periodic (OL) potential 
with strength Vb and period d — i^jk. In the following analysis below, we adopt the following representative values for 
parameters of the potential: = 0.21, = 0.5 and k = 0.3, in which case the four smallest eigenvalues of operator 
L are found to be 

cJo = 0.3585, cji = CJ2 = 0.3658, c^s = 0.3731. (4) 

In a weakly nonlinear setting, we implement the natural possibility of a four-mode approximation, based on a 
Galer kin- type expansion of u{x^y^t) and truncation of the higher-order modes. We denote the ground and first three 
excited eigenstates of the linear operator L, shown in Fig. (TJ as uq and i/i,2,3- This set constitutes a natural minimum 
basis for the Galerkin approximation in the system of four potential wells coupled by tunneling. Eigenstates uj 
(j = 0, 1, 2, 3) can be chosen to be real, given the Hermitian nature of the operator L. Then, solutions of Eq. ([1]) for 
values of the chemical potential in a vicinity of linear eigenvalues may be approximated by linear combinations 
of the four eigenfuctions. 

Actually, it is more convenient to use a transformed basis, {<I>o, ^2, *^3}, as shown in Fig. [21 which is based on 
populations of the four wells. This basis is generated from the original set by a linear transformation: 

[^0 ^2 ^3] = [^0 ui U2 Us] T, (5) 

where the appropriate transformation matrix is 



(6) 



Each mode {j = 0, 1, 2, 3) is localized in one of the four wells, with the four of them constituting an orthonormal 
set. 

Using the new basis, we can readily reformulate the four-mode decomposition as 

3 

u{x,y,t) = ^Cj{t)^j{x,y), (7) 
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FIG. 1: (Color online) The wave functions of the ground state, no, and the first three excited states, ui, U2 and us, for the 
four- well potential of Eq. (|3]) with Q = 0.21, Vb = 0.5 and k = 0.3. Note the difference in the grayscale (color, in the online 
version) bars in the first and three others panels, related to the fact that the wave function of the ground state is positive, 
while the excited states feature sign-changing patterns. 
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FIG. 2: (Color online) Basis modes {$o, ^i, ^2, ^3} that are localized in each of the wells. 



with time-dependent complex amplitude Cj(t), j = 0,1,2,3. Substituting Eq. ([7j) into Eq. ([T]) and projecting onto 
the orthonormal basis {<l>o, ^3}, we derive, by means of straightforward algebra, the following system of four 

ordinary differential equations (ODEs), 

icj = ujj + sAj\cj\^Cj -^s'^Bjk {2\ck\^Cj + c^c*) 

+ s ^[Dkj\ck\^Ck + Djk {2\cj\^Ck + c^jcl)] + s ^ Ekji (2|cfcpQ + 4^*) 

ky^j k^l^j^k (8) 

+ 5 ^ Ejkl {c*CkCi + CjClci + CjCkC*i) + s G ^ 

^^^^^^^ k^i^m^k 

with the summation performed over k,l,m = 0,1,2,3. To cast these equations in a more compact form, we have 
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(9) 



where 70 = c^o + cc;i + CJ2 + ^3, 7i = c^o + — CJ2 — ^3, 72 = cc;o - cc;i + CJ2 - ^3, and 73 = uq — ui — U2 ^ 
Notice tiiat, for tiie underlying eigenvalues (|4|), one has 71 = 72 = ^0 — ^3 and 73 0. Furthermore, Eqs. (|8|) involve 
nonlinear coefficients given by overlap integrals, viz.^ An = jj^^dxdy^ Bmn = JJ^m^n^^^V^ Dmn = JJ^m^ndxdy, 
Eimn = JJ^f^m^ndxdy, G = JJ 1^2^ 3 dxdy ^ with l^m^n = 0,1,2,3; these indices must be mutually different 



wherever they appear in the coefficients before the nonlinear terms. 

For our choice of the parameters of the potential, the overlapping between modes is weak (see Fig. [2]), therefore 
all the overlap integrals are much smaller than the A^s. Neglecting these small overlap terms leads to the following 
simplification of Eq. ([8]): 



■Aj\cj\ Cj, j 



0,1,2,3. 



(10) 



It has been checked that the latter reduction of the four-mode equations very slightly affects the accuracy of the 
solutions, while it renders the identification of various bifurcation branches significantly easier. Furthermore, this 
reduction is more convenient in simulations as we may use these solutions as inputs for generating numerical solutions 
of the full GP system, as explained below. 

In this 2D setting, we seek both real and complex stationary solutions to the ODE system. Substituting Cj{t) = 
pj{t)e^'^^^'^^ into Eq. ([10]), we split them into real equations for pj and Lpj\ 

(11) 

)], (12) 
(13) 

)], (14) 

with the equations for p2,3 and (/P2,3 obtained by interchanging the indices, < — > 2 and 1 < — > 3, except for in 70 and 
71- 

Looking for solutions with constant pj and ipj which are integer multiples of tt, we reduce Eqs. (fTT]) - (p!4|) to 
a set of four algebraic equations for pj, which can be used to derive a complete set of stationary solutions of the 
four-mode truncation. These were further used as initial guesses to find numerical solutions of the full system of the 
GP equations. Moreover, our analysis of the four-mode system indicates that nontrivial complex solutions in this 
setting are only possible in the form of discrete vortices^ i.e., solutions with phase sets ipj = nj j = 0, 1,2,3 [22[, 
which have been studied in detail in Refs. [19] and [21] (we also briefly consider them here). 



Po 


= ^71 [Pi ^H^i -^0) H 


- ps sm{(ps - (po)' 


7 








^0 


= {P- ^7o) - sAopI - 


1 rPl / 

-7i[— cos((/?i - 
4 Po 


^0) - 


Po 




- ^0 


pl 


= ^7i[Posin((/?o - 9^1) H 


- P2 sm{(p2 - (pi)\ 










^1 


= {P- ^70) - sAipl - 


1 .Po ( 
-7i[— cos((/?o - 
4 p\ 


^1) - 


pl 


COs((/92 


- ^1 



III. NUMERICAL RESULTS 



A. Attractive interactions 



Let us first present results of numerical simulations pertaining to attractive interactions (alias self-focusing nonlin- 
earity) case, i.e., s = —1 in Eq.([T]). Our basic bifurcation diagram, shown in Fig. [3l displays the squared norm of the 
solution (which physically describes the number of atoms in BECs or the power in optics), = JJ\u{x, y, t)p dxdy, as 
a function of the chemical potential p. The left panel of Fig. [3] presents the full numerical bifurcation diagram, which 
involves twelve real and one complex solutions (for the latter solution, N is the same as one of the real branches, hence 
this branch is not visible as a separate curve in the diagram). The companion diagram in the right panel is obtained 
from the above-mentioned algebraic system for stationary solutions produced by the the four-mode reduction, and 
demonstrates good agreement with its numerical counterpart. 

The twelve real branches are labeled mainly according to their relation to the populations of the four wells. To 
support our explanation, we introduce a symbolic representation that we developed in the form of 2 x 2 matrices. 
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FIG. 3: (Color online) Top panels: squared norm N (normalized number of atoms in BECs or power in optics) of numerically 
found solutions of Eq. ([T]) (left), and their counterparts predicted by the four- mode approximation (right), for attractive 
interatomic interactions (s = — 1), as a function of /x, i.e., respectively, the chemical potential or propagation constant. The 
bottom panels are segments of the top left panel. The (blue) soUd lines and (red) dashed lines denote stable and unstable 
solutions, respectively. The branches are explained in the text and their profiles and stability are detailed in Figs. |4][7l 



labeling different waveforms that arise in the diagram, as follows: Al 

C4: = ^ ^ and Dl = ^ ^ ^ ^ . In this representation, 1, —1 and have the obvious meaning, 

by indicating that a particular well is or is not populated, and the phase of the wave function in it being or tt in the 
cases of +1 and —1, respectively, when populated. Symbol e denotes either a small (but nonzero) population in one 
of the wells, or a symmetry-breaking effect (when some of the density peaks feature values ±1 ±6, thus being slightly 
different from ±1). The labeling is then defined as follows: branches A1-A4 have the same amplitude at the four wells 
as long as they are populated, branches B1-B3 feature two pairs of peaks with different amplitudes, branches C1-C4 
have three different amplitudes, while Dl has all of its four peaks different. The waveforms in the top rows of Figs. 
HJIIl display prototypical realizations of the relevant branches. Their stability properties are illustrated, as a function 
of eigenvalue parameter /i, in the bottom rows. 

We will now explain in detail solutions appearing in the full bifurcation diagram, starting from the linear limits 
(A^ ^ 0). First, we look at the group of solutions related to branch Al, as shown in the bottom left panel of Fig. [3l 
This branch arises from the symmetric linear mode at /i = cjq, i.e., the ground state in the linear limit, uq. Accordingly, 
Al features four identically populated wells. The analysis demonstrates that it is stable near the linear limit, but 
is soon destabilized, due to the emergence of branches CI and Bl, through subcritical and supercritical pitchfork 
bifurcations, respectively, around /i = 0.355. In other words, there are two consecutive steady-state bifurcations, 
in the language of Ref. [29], in two different subspaces, in which one unstable solution, CI, collides with Al and, 
simultaneously, a pair of eigenvalues emerges on the real axis for the Al branch with the decrease of /i (in a subcritical 
pitchfork); then, a super-critical pitchfork takes place in another subspace, in which Bl retains only one real pair, 
while another pair passes through the origin along the Al branch. The actual "pitchfork" cannot be visualized here 
in the usual manner, because any of the four equivalent versions of Bl (obtained by 7r/2 rotation) have the same 
value of A/", and so are represented by the same branch in the graph. Branch Bl, which is unstable due to a pair of 
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FIG. 4: (Color online) Top panels show the profiles of wave functions of branches Al, A2, A3 and A4 (from left to right) 
at /i = 0.34. The bottom panels display real parts of the unstable eigenvalues of the respective branches as functions of the 
parameter /x. 
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FIG. 5: (Color online) Top: profiles of wave functions that represent branches Bl, B2 and B3 (from left to right) at fi — 0.34. 
Bottom: real parts of unstable eigenvalues of the corresponding branches as a function of /i. The dashed-dotted line in the 
bottom right panel indicates a complex quartet of eigenvalues. 



real eigenvalues in the linearization around it throughout its domain of existence, features two of the wells on one 
side being less populated than the other two. Configuration Bl becomes increasingly more asymmetric as it deviates 
from Al. A noteworthy feature is shown by branch CI, which bifurcates from Al at almost the same place as Bl: 
after having emerged, it tends to be located on the left of Al as are all other branches bifurcating from Al. However, 
within a narrow interval of /i, its norm decreases {dN/dfi > 0) slightly before starting to grow as usual {dN/dji < 0). 
Naturally, when the norm decreases the solution is destabilized and then it remains stable after the turning point, 
which is explained by the well-known Vakhitov-Kolokolov criterion [2^. Branch Al is endowed with two identical 
pairs of real eigenvalues by Bl and CI upon their bifurcation (which is shown as the dashed line in the bottom left 
panel of Fig. [H) As N grows further, a subsequent bifurcation, at /i = 0.3519, leading to the emergence of branch 
B2, adds yet another real eigenvalue pair to Al; this means that Al possesses in total three real eigenvalue pairs for 
sufficiently large vaues of N. Branch B2 features two wells on the diagonal which are less populated than the other 
two, and it is unstable, with two pairs of real eigenvalues, near the point where it is generated by the bifurcation 
from Al; however, one of the pairs is eliminated by the emergence of a new branch, C2, from B2 through a pitchfork 
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shortly afterwards. Branch B2 then remains unstable with one real pair, while C2 (with three principal sites, one of 
which is of lesser amplitude than the other two) is unstable due to two real eigenvalue pairs. 

This description encompasses all branches of stationary solutions which can be traced back to the ground state of 
the linear system. Detailed information for the wave function profiles and the development of the real eigenvalues 
associated to them is presented in Figs. |4]l6l 






FIG. 6: (Color online) Top: profiles of wave functions of branches CI, C2, C3 and C4 (from left to right) at /i = 0.34. Bottom: 
real parts of unstable eigenvalues of the corresponding branches as a function of the eigenvalue parameter /i. 




FIG. 7: (Color online) Left: the profile of wave function of branch Dl at /i = 0.34. Right: real parts of its unstable eigenvalues as 
a function of the eigenvalue parameter /i. The dashed-dotted line in the right panel indicates a complex quartet of eigenvalues. 

Next we turn to the states originating from the second linear mode, as shown in the bottom right panel of Fig. [3l 
Branch A2 starts from the respective eigenvalue, fi = uji = cj2, which pertains to the first and second (degenerate 
at the linear limit) excited states. This branch emerges as an unstable one, carrying a real eigenvalue pair. The 
respective wave function profile features four wells populated with the same amplitude but tt out-of-phase between 
the two sides, see Fig. HI Branch B3 emerges from A2 through a supercritical pitchfork at /i = 0.3623, lending 
another real eigenvalue pair to A2. Similar to the case of Bl (as it separates from Al), in the B3 state two wells on 
the one side tend to be less populated than the other two, as this branch moves further from A2. Branch B3 remains 
unstable with one real eigenvalue pair, until getting stabilized by another pitchfork bifurcation which takes place at 
/i = 0.3589; this simultaneously gives rise to a new unstable branch, Dl, that has different populations in all four 
wells. Notice that both B3 and Dl pass through a Hamiltonian-Hopf bifurcation (alias 1:1 resonance, in terms of Ref. 
[29]), which means that, in the relevant parametric interval (0.3362 < ji < 0.348 for B3; 0.3444 < ji < 0.3553 for Dl), 
B3 is destabilized by a complex quartet of small eigenvalues in the linearization around the stationary solution, while 
Dl remains unstable, but with one real eigenvalue pair and a complex quartet (the dashed-dotted lines in the bottom 
right panels of Figs. [5] and [71 refer to this effect). 

Furthermore, branch A4 bifurcates from the same linear mode as A2. It is unstable near the linear limit due to a 
Hamiltonian-Hopf bifurcation, but with the increase of N it becomes stable. Branch A4 features a waveform in which 
only two wells lying on the diagonal are populated, with the same amplitude but opposite signs. 
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Branch A3 arises from the third excited hnear mode at /i = CJ3. In this case, two wehs on the diagonal are populated 
with equal amplitudes, while in the other two the amplitudes are of opposite signs. It is the unique stationary solution 
which remains stable across the entire bifurcation diagram (despite the fact that it has three pairs of purely imaginary 
eigenvalues with negative Krein signature [sO*], which in principle, can give rise to Hamiltonian-Hopf bifurcations). 

Finally, branches C3 and C4, which are located slightly below A2 in Fig. [3l correspond to a pair of states which arise 
through a saddle-node bifurcation at some critical value of chemical potential (/i ^ 0.348, for the model's parameters 
chosen in the present case.) Branch C4 (the one with higher values of N) is unstable with a real eigenvalue pair, while 
C3 remains stable, except inside a short instability interval, which is accounted for by a Hamiltonian-Hopf bifurcation. 

In addition to the above real stationary states, we have also found complex solutions in the form of vortices l^,^. 
A typical example of such a solution is shown in Fig. [8l Throughout the regime of parameters considered herein, 
such solutions have been found to be linearly stable. 




FIG. 8: (Color online) The absolute value (left) and phase (right) of a vortex state for /i = 0.34. 



It is interesting to note that, for all the solutions considered herein, in the large- limit their stability characteristics 
coincide with what can be suggested by the DNLS model considered in Ref. [22] (see also Refs. [24, 25] for the 
corresponding ID and 3D stability results). Gross features of these findings are that, whenever two adjacent sites are 
in-phase, a real eigenvalue pair is expected to emerge due to their interaction, while whenever such sites are out-of- 
phase, the relevant eigenvalue is expected to be imaginary [26], but with negative Krein signature \22\, which implies 
a potential for a Hamiltonian-Hopf bifurcation. It should also be noted that, in the limit of the infinite lattice, it is 
naturally expected that the asymmetries observed herein in many of the branches will disappear (i.e., the amplitudes 
in different wells will be equal) - see also a relevant discussion in Ref. [l^ • 



B. Repulsive interactions 



We now briefiy discuss the case of repulsive interactions (alias self-defocusing nonlinearity ) , corresponding to 5 = +1 
in Eq. ([1]), with an objective to highlight its similarities with and differences from the case of attractive interactions. 
The bifurcation diagram for the model is displayed in Fig. [9l 

The solutions are labeled so as to match the self-focusing case, by means of the appropriate staggering transformation 
[2^. The latter effectively converts the defocusing nonlinearity into a focusing one by changing the relative phase of 
nearest-neighbors from to tt and vice versa, while preserving the relative phase of next-nearest-neighbors. In this 
way, each solution in the defocusing case is linked to its counterpart in the focusing model through this transformation. 
Following this relation, and adopting the same matrix symbolic representation used for the focusing case in section 

nil A( the branches of solutions are labeled as follows: Al = ^ ^ ,^42=^ | 1^ ' ^ (^1 ' ^ 0^ 1 

Thus, in this case, the symmetric ground state of the system is A3, which is stable for arbitrary values of N. Branch 
A2 is immediately unstable, starting from the linear limit. B3 bifurcates from A2 and remains unstable before getting 
stabilized through giving birth to Dl (and then becoming destabilized again). Branch Al is stable near the linear 
limit, but is subsequently destabilized due to bifurcations that give rise to Bl and CI, and an additional real eigenvalue 
pair arises at higher value of N due to the emergence of B2, from which another new branch, namely C2, arises in 
turn. Branches C3 and C4 exist for a while (when N is large enough), and then collide at /i = 0.389 (for the values 
of parameters adopted herein). The types of the bifurcations, the emergence of the corresponding solutions, and the 
corresponding stability properties were found to be in direct correspondence to the case of self-focusing nonlinearity. 
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FIG. 9: (Color online) The norm of the numerical (left) and approximate (right) stationary solutions to Eq. ([T]) with the self- 
defocusing nonlinearity (s = +1), as a function of the parameter /i (the chemical potential in BECs or propagation constant in 
optics). The labels of the branches are explained in detail in the text. 

provided that one takes into account the staggering transformation relating the repulsive and attractive case models 
as indicated above. 



C. Dynamics 

We now proceed to investigate the evolution of unstable states in the model with the self-focusing nonlinearity. 
To this end, for each unstable branch, a small perturbation is added to the most unstable eigendirection of the 
linearization near the original stationary solution, at /i = 0.335. Results of the simulations are presented in Fig. [TOl 

Panel (a) shows the behavior of solution Al, which, as a result of the instability, starts oscillating between a state 
where all four wells are populated and one in which only two diagonal wells are not empty. Panel (b) depicts a 
periodic oscillatory behavior of A2, also between four and two populated sites, but in this case the continuously 
populated sites are adjacent to each other. Unstable mode A4 (panel (c)) features only two non-empty wells, with the 
symmetry-breaking instability resulting in the enhanced population of one of the two. As explained in section UlI A[ 
modes Bl, B2 and B3 have two very weakly populated wells, in comparison with the other two. Since we employ 
isosurface \u {x^y^t) = k (where k is the half of the maximum density at t = 0) to plot the dynamics in Fig. [TOl the 
evolution in the weakly populated wells is not visible (which indicates that they play a minor role in the dynamics). 
Panel (d) shows that mode Bl sustains a symmetry breaking similar to that of A4, but between adjacent sites. On 
the other hand, modes B2 and B3 appear to be oscillating between the two dominant wells roughly periodically, as 
shown in panels (e) and (f). Mode C2 [in panel (g)] oscillates between three and two populated sites (the seemingly 
empty well is actually a weakly populated one, similarly to the modes of type B, see above). Mode C3, whose weak 
instability is caused by a quartet of eigenvalues, is also "breathing" within the respective set of three predominantly 
populated wells, as shown in panel (h). Finally, mode C4 [shown in panel (i)] involves a complex symmetry-breaking 
pattern, with different numbers of wells populated at different times, while mode Dl (panel (j)) oscillates between 
three- and two- well asymmetric configurations. 

IV. CONCLUSION 

In this work, we have studied stationary and dynamical properties of the two-dimensional nonlinear 
Schrodinger/Gross-Pitaevskii equation, including a four- well external potential, with both signs of the nonlinearity, 
self-attractive (focusing) and self-repulsive (defocusing). The model applies to a BEG confined in a highly anisotropic 
harmonic trap, with the transverse confining frequency being much smaller than the axial one, which results in a 
planar configuration. In this context, the four- well potential can be generated by a combination of the transverse part 
of the harmonic trap and an optical lattice. The same model may describe the propagation of an optical beam in a 
bulk nonlinear medium with an embedded four-channel guiding structure. 

In our analysis, first we developed a four- mode approximation, which strongly simplifies the identification of sta- 
tionary solutions. Using this approximation, we were able to find the four (two of which are identical) symmetric 
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FIG. 10: (Color online) The spatiotemporal evolution of unstable states, represented by the respective density isosurface, 
2/, = /c, where constant k is taken as half the maximum value of the density distribution at t = 0. The results are 
arranged as follows. Top panels: Al, A2, A4; middle panels: Bl, B2, B3; bottom panels: C2, C3, C4, Dl. 

and antisymmetric linear modes, and all branches of asymmetric solutions emerging from them in the model with 
the self-focusing nonlinearity. The linear-stability analysis demonstrated how pitchfork and saddle-node bifurcations 
change the stability of the branches. We have shown that in the limit of strong nonlinearity, properties of localized 
modes in the model with either sign of the nonlinearity can be, roughly, understood on the basis of earlier results 
pertaining to the corresponding discrete NLS model. We have also described the evolution of all unstable solutions, 
observing, typically, the emergence of symmetry-breaking instabilities and the emergence of respective oscillating 
solutions. 

It would be interesting to investigate how these four-site configurations may be embedded into a larger potential 
pattern, with 9 or 16 wells, and examine whether the symmetry-breaking bifurcations considered above are sustained 
(or how they are modified) within the larger pattern. In this context, a conjecture that would be worthwhile proving 
is that, in an infinite periodic lattice formed by potential wells, the nonlinearity can support 2D solitons and localized 
vortices with various symmetries, but not confined asymmetric states. This conjecture is suggested by results reported 
for infinite linear [31] and nonlinear [32] potential lattices. 
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